Two technologies are widely used and accepted as reliable methods for performing bathymetric surveys. These technologies are considered the gold standards of bathymetric surveying. One technology is ship-based side-scan sonar, and the other is airborne bathymetric LIDAR (light detection and ranging). Both meet the industry IHO 1 standard in depth accuracy and resolution. They are however expensive and time consuming. There is a need for less costly surveys with more rapid turnaround.
Two cost-saving alternative bathymetry surveys technologies have received the most attention and development efforts. The first is photogrammetric bathymetry, basing the depth on the relationship to the radiance of the ocean bottom. The main problem with the photogrammetric method is the need to know the optical properties of the water and ocean bottom. Usually these are not well known and can be variable over a small area. The second alternative technology is “wave-kinematics bathymetry” (WKB), for which systems and algorithms are discussed here. However, despite decades of research and development, neither technology has achieved much acceptance and commercial success.
WKB systems and algorithms use the relationship between the velocity of gravity waves and depth. In the simplest terms, the speed of ocean waves decreases as waves propagate into shallow depth and “feel” the bottom. More precisely, the frequency-wavenumber (ω,k) spectrum transforms according to the well known linear gravity wave dispersion formula,ω0({right arrow over (k)}|d,{right arrow over (u)})=√{square root over (g|k|tan h(|k|d))}−{right arrow over (k)}·{right arrow over (u)}, where d is the ocean depth and {right arrow over (u)} is the surface current vector. This equation is the gravity wave dispersion approximation for wave heights, H, that are small relative to the depth. While the linear gravity wave approximation is the one most often used, alternative approximations have been proposed for moderate wave heights, H→d. An example is the empirical approximation [Hedges, 1976]:ω0({right arrow over (k)}|d,{right arrow over (u)})=√{square root over (g|k|tan h(|k|d(1+αH|d)))}−{right arrow over (k)}·{right arrow over (u)}, where α is an empirically derived constant in the range of 0.4-0.5 [Flampouris et al., 2009] and H is the wave height (either known or guessed-typically ˜1 m). It can be readily seen that this is the linear gravity wave dispersion with a bias correction. Note that WKB algorithms are useful primarily for d>H. The region where d≦H defines the surf zone where waves break and become bore waves. This region is typically excluded from bathymetry data. As such, the non-linearities that occur as H→d are of little practical significance to real-world bathymetry and are included herein only for completeness.
Several WKB sensor systems have been implemented over the past 25 years. They are mostly experimental, but a handful offer bathymetry surveying as a commercial service. Various imaging systems can be used.1 In all such systems, a sequence of images is recorded to observe the wave motion over a dwell time, T. The image data is thus referred to as an image cube comprised of three dimensions, two spatial dimensions and time. 1 The purest image data for bathymetry would be images where the pixel intensities are linearly proportional to the ocean surface heights. Most sensor systems do not actually image surface height. The most common sensors used in this method are visible light cameras and microwave radars. Cameras image the ocean waves by reflection of sun and sky light; microwave radars image waves through Bragg scattering from wind generated capillary waves. Other sensor modalities exist for imaging waves. Some examples: thermal IR, as reported in [Farber, 1995] and [Dugan, 1996], and topographic LIDAR [Wang, 2002]. The topographic LIDAR (not to be confused with bathymetric LIDAR) is the only sensor providing direct measure of the surface height. All the others do not but under appropriate conditions can be reasonable approximations of height sensing.
Two types of algorithms used to process the image cubes into bathymetry have been previously described:                3D Fourier Analysis: The image data is partitioned into “tile-cubes” of XY in space, and T in time. The cubes are transformed from three-dimensional space-time data into three-dimensional wavenumber-frequency spectra. These algorithms are the most common and most practiced in the literature. This technique was first mentioned in [Hoogeboom, 1986] and refined in many subsequent works, e.g., [Dugan et al., 2000]. 3D Fourier Analysis uses very large image cubes.        2D Fourier Analysis: The “tile-cubes” are Fourier transformed into two-dimensional wavenumber spectra. The advantage of this approach relative to 3D Fourier Analysis is that the dwell time can be shorter. The method, including some experimental results, are described in [Abileah, 2006A], [Abileah, 2006B], and [Abileah, 2010].        
The primary output of these methods is typically a bathymetry chart, sometimes also referred to as a nautical depth chart. It is noteworthy that all WKB methods can also provide a map of the ocean surface current. For bathymetry, the current is considered a nuisance parameter, not the primary product of interest. However, some users are actually more interested in mapping currents, for example to detect rip tides or river outflows ([Abileah, 2007]).
For the end users of any bathymetric survey, the important considerations are cost, turnaround time, depth accuracy, and spatial resolution. With respect to the first two criteria, cost and turnaround, the known WKB systems compete very well with other bathymetric methods. However, depth accuracy and spatial resolution are found wanting in the 2D and 3D Fourier Analysis methods. The resolution is so poor that even with the significant cost benefit, WKB has attracted little commercial interest. Other, much more costly technologies (side scan sonar, airborne bathymetric LIDARs) are used for most actual commercial bathymetry surveying.
The reason for the poor performance has to do with a fundamental limitation of Fourier transforms. An excellent discussion of the issues can be found in [Piotrowski et al., 2002]. Briefly, the problem is that the depth accuracy with these methods is proportional to the resolution of the wavenumber-frequency transform, which in turn is inversely proportional to the image cube dimensions, the length on a side of the XY tiles, and the dwell time, T. The 2D and 3D processing algorithms are a compromise between depth accuracy and resolution.
For example, suppose that 3D Fourier Analysis is performed using tile-cube dimensions of 250 m×250 m×100 s and achieves a respectable 10% accuracy in depth. If the spatial dimensions are reduced to 25 m×25 m to get better spatial resolution, then according to Piotrowski et al., the depth error blows up by a factor 10 (from 1/250 to 1/25) and the error is then 100% of depth. A bathymetry chart with 100% error is not useful. 2D Fourier Analysis can achieve good depth estimates with shorter dwell time, and the spatial resolution can be slightly better than for 3D Fourier Analysis, but the depth accuracy is still limited by the XY tile dimensions used for processing.
In summary, the mathematical approach used in 2D and 3D Fourier Analysis cannot achieve high depth accuracy and spatial resolution simultaneously. For this reason, WKB has had limited commercial success up to now. However, it is still useful to understand how the methods work, as the instant invention uses related method steps and concepts.
3D Fourier Analysis: Several versions have been described; all use a three-dimensional Fourier transform, but other details may differ. The following description is closest to the implementations in [Hoogeboom 1986], [Dugan et al., 2000], and [Trizna, 2001].
A sequence of raw images is obtained. Pre-processing removes image sensor artifacts (e.g., bad pixels) and orthorectifies the images into fixed-Earth Cartesian coordinates using standard methods. Next, the images are partitioned into tiles of dimensions X by Y, usually square. The spatial dimension is driven, as mentioned above, by the resolution vs. accuracy tradeoff. For purpose of illustration, assume that the tiles are 250 m×250 m and that the full image area is 5 km×5 km. There are thus (5000/250)×(5000/250)=400 tiles. If the tiles are overlapped 50% in each dimension (meaning that tile centers are separated by 125 m in each dimension), there will be (5000/125)×(5000/125)=1600 tiles. A calculation is then made for each tile to obtain a best fit to the gravity wave dispersion formula using average depth of the tile as a fitting parameter. The calculation produces one depth solution per tile. That is, the entire XY area is represented with one depth value; there will be 400 depth estimates over the entire 5 km×5 km image with non-overlapping tiles or 1600 depth estimates with overlapping tiles.
Let S represent the space-time cube of N images {S1({right arrow over (x)}), S2({right arrow over (x)}) . . . SN({right arrow over (x)})} for one tile. The images are taken at intervals 2 so that the time span is T=(N−1)τ. The cube is transformed into a wavenumber-frequency power spectrum with a 3D Fourier transform, ℑ3:P({right arrow over (k)},ω)=[{S}]|2.
In the transform space, wave energy is found to be concentrated in a thin wave dispersion shell which can be easily separated from non-wave energy [Nieto-Borge and Soares, 2000]. The depth, d, and the ocean current vector, {right arrow over (u)} are determined by finding the best fit of the gravity wave dispersion equation, ω0({right arrow over (k)}|d,{right arrow over (u)}), to the measured power spectrum, P.
The best fit is usually found by seeking to minimize/maximize some objective function J on the power spectrum. For example, [Tang et al., 2008] suggested using
            argmin              d        ,                  u          →                      ⁢          J      Tang        =            ∑                        k          x                ,                  k          y                      ⁢                  W        ⁡                  (                      k            →                    )                    ⁢                        ∑          ω                ⁢                                            (                              ω                ±                                                      ω                    0                                    ⁡                                      (                                                                                            k                          →                                                |                        d                                            ,                                              u                        →                                                              )                                                              )                        2                    ⁢                                    P              ⁡                              (                                                      k                    →                                    ,                  ω                                )                                      .                              
W({right arrow over (k)}) is not absolutely essential, but is recommended to improve the accuracy of the solution. An appropriate W can increase the signal-to-noise ratio (SNR) of waves by assigning greater weight to the wavenumbers with ocean wave energy and less weight where there is more noise than waves based on knowledge of where such energy is likely to be present. A separate algorithm is required to derive the optimum W (see [Abileah, 2010] and [Tang et al., 2008]).
2D Fourier Analysis: As in 3D Fourier Analysis, a sequence of images are pre-processed through ortho-rectifications to fixed-Earth Cartesian coordinates, and the image space is partitioned into XY tiles. Processing proceeds one tile at a time. The 2D Fourier transform ℑ2 is applied to each of the N images, producing the 2D complex spectra,Fn({right arrow over (k)})=[Sn({right arrow over (x)})].  (1)
A wave propagation kernel is defined as,Φ±(kx,ky|d,{right arrow over (u)})≡e±iω0(kx,ky|d,{right arrow over (u)})τThe sign in the exponent depends on whether the waves are traveling in +x or −x. Assume that the xy coordinates are rotated such that all waves propagate approximately parallel to the x-axis. Whether the primary wave direction is +x or −x is usually obvious by inspection: waves normally travel towards the shore. To automate the algorithm one can also determine the correct sign based on a test of the inter-frame coherence. The condition(Φ+Fn)*Fn+1>(Φ−Fn)*Fn+1  (3)implies that waves are propagating in +x and Φ+ should be used. If the direction of the inequality is reversed, use Φ−. In the open ocean one often finds different waves systems (e.g., swell waves and wind waves) traveling in opposite directions. This is unlikely to occur near shores where all waves travel more or less towards shore. For purpose of bathymetry, the complication of opposing wave systems is not important.
For gravity waves, the Fourier transform at time step n is simply related to another at time step n+m by the equationFn+m=Φ±mFn.
Using this relationship, depth can be determined from a pair of images, {S1({right arrow over (x)}), S2({right arrow over (x)})} by minimizing the objective function
            argmin              d        ,                  u          →                      ⁢          J      Δ        =            ∑                        k          x                ,                  k          y                      ⁢                  W        ⁡                  (                      k            →                    )                    ⁢                                                                              F                1                            -                                                Φ                                      -                    m                                                  ⁢                                  F                  2                                                                          2                .            
This algorithm was first described in [Abileah, 2007]). Experimental results on real data were reported earlier in [Abileah, 2006A, 2006B].
The extension to more than two images is simply a sum over a multiple of image pairs ([Abileah, 2010]),
            argmin              d        ,                  u          →                      ⁢          J      Δ        =            ∑                        k          x                ,                  k          y                      ⁢                  W        ⁡                  (                      k            →                    )                    ⁢                        ∑                      n            =            1                                N            -            m                          ⁢                                                                        F                n                            -                                                Φ                                      -                    m                                                  ⁢                                  F                                      n                    +                    m                                                                                            2                    History
The physics of how ocean waves are transformed in shallow depth has been known for over a century. The earliest literature suggesting the application to bathymetry, [Biesel, 1952], [Hart and Miskin, 1945], [Williams, 1946], and Fuchs, 1953], was written long before adequate computing and sensor technology existed. The key enabling technologies for wide use of WKB came in the 1980s. These were: the fast Fourier transform (invented in the mid 1960s), an appreciation of the use of Fourier transforms for analysis of ocean waves (about 1980), and increasingly more effective digital imaging systems (1980 to present).
The Fourier transform is needed to compute the frequency-wavenumber spectrum. Theory predicts that energy of linear gravity waves concentrates in frequency-wavenumber space on a thin shell or manifold. Two seminal papers appearing in the 1980s, one on visible light imaging [Gotwols and Irani, 1980] and the other on marine radar imaging [Young et al., 1985], proved the theory correct in relatively deep water.
Right on the heels of [Young et al.], a Dutch team tested a marine radar on the dunes of Scheveneingen beach, imaging shoaling waves. With their initial data they demonstrated that the wave dispersion spectrum shifted from deep water to shallow as predicted by theory [van Halsema and Kleyweg, 1986 & Hoogeboom, 1986]. The Hoogeboom paper at the IGARSS'86 conference in Zurich was the first public presentation of WKB with the 3D method.
Others repeated and improved on the Dutch success in the following 25 years (from 1986 to the present), all using some variation of the 3D algorithm. Table 1 summarizes the major research efforts. The table also lists the tile-cube dimensions used. The cube dimensions are not always explicitly stated in papers and in some cases have to be deduced from several sources. Some of the bathymetry results by the GKSS (Germany) researchers are based on their so called DiSC algorithm that does not tile and is claimed to produce better spatial resolution than achieved with the traditional tile method.2 2 The DiSC algorithm is claimed to improve resolution over the conventional 3D tiles. It starts with a 3D Fourier transform as in other 3D methods, then inverse Fourier transforms a narrow frequency range into XY [Sent et al., 2008]. The idea is that a single frequency wave shifts to higher wavenumber as it propagates into shallower water. Some limited DiSC test results are presented in the literature and show neither improved resolution or depth accuracy.
TABLE 1History of 3D WKBPrincipalInvestigator/Tile-cube dimensionsOrganization/SpatialTemporalSensorYearReference(meters)(seconds)Hoogeboom/TNO,1985van Halsema and 750100the Netherlands/Kleyweg, 1986Marine X-bandHoogeboom, 1986RadarSmith/DSIR, New 1991Smith et al., 199132060Zealand/DopplerRadarDugan/Arete/1995Dugan et al., 1996~100070Airborne IRFarber et al., 1995cameraBell/Proudman1999Bell, 1999~300160OceanographicBell et al., 2001Laboratory, UK/Bell et al., 2004Marine X-bandradarHessner/GKSS/1999Hessner et al., 1999300 × 600121WaMos II radarDugan/Arete/2000Dugan et al., 2000256128AROSS airbornePiotrowski and cameraDugan, 2002Dugan et al., 2001GKSS, Germany/2003Dankert, 2003N/A128WaMoS radarGKSS Germany/2005-Flampouris et al., 2009N/A460WaMoS II radar2009Flampouris, 2006Senet et al., 2008]GKSS Germany/2009Hessner and Bell, 2009240600WaMoS II radar
Another approach to WKB has been developed specifically for conditions where the predominant waves are swell. The unique features of swell is that the waves are highly coherent, and the structure of crests and troughs can be easily identified and tracked over a time period, either by eye, or by some standard image processing algorithm such as optical flow. Theory tells us that the wave frequency is constant for a wave propagating from deep water to shore. Only the wavelength, and velocity (wavelength and velocity being related by a simple formula) change with depth. The wavelength can be measured by the distance from one crest to the next, and velocity can be measured by the displacement of a wave crest in the time interval between two images. Very precise depth estimates and good resolution is possible. There are no Fourier transforms and the aforementioned problems with Fourier transforms are eliminated. One example of this principle put into practice is the 1985 work of scientists at the Naval Research Laboratory-Stennis Center with a pair of images taken six seconds apart with an aircraft reconnaissance camera [Caruthers et al., 1985]. The wave crests in each image were manually traced and recorded and depth determined from the apparent crest displacement in the six-second interval. In actual fact, the wave conditions in the Caruthers et al. data were not ideal swell waves (rather, they were a combination of swell and wind waves), which partially explains their 30-40% depth errors. (The other explanation is errors from the manual tracing of wave crests.) More recent work with this technique is given by [Stockdon and Holman, 2000] and [Takewaka, 2000] using full digital imaging and data processing. All these works are unable to separate the effect of depth and currents. They implicitly assume that the current is zero in the direction of wave propagation. If the current is not zero, the estimated depth is biased.
The 2D approach was motivated by the introduction of the IKONOS satellite, launched in 1999. Before the IKONOS satellite, commercial satellites were generally taking only single images of an area of interest. IKONOS had a camera re-pointing capability such that the same area could be imaged 6-7 times during an overpass, at intervals of approximately 13 s. It was then possible to consider using IKONOS imagery for WKB, but the conventional 3D algorithms could not be used. The imaging geometry changes quickly in satellite orbits, limiting the effective dwell time to the interval between two images (13 s). For such short dwell times the rms error of the 3D method is 100% of depth.
The 2D approach solved this problem, and an early version of the 2D process was applied to pairs of IKONOS images of the San Diego inlet area [Abileah, 2006A, 2006B, & 2007]. An updated version of the 2D was applied to X-band radar images obtained by Dr. Dennis Trizna with his experimental radar at the US Army Corp of Engineers Field Research Facility (FRF) [Abileah, 2010].
It is apparent, however, that the mathematical approaches used in 2D and 3D Fourier Analysis cannot achieve high depth accuracy and spatial resolution simultaneously.